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Abstract 

Rich and complex time-series data, such as those generated from engineer- 
ing systems, financial markets, videos or neural recordings, are now a common 
feature of modern data analysis. Explaining the phenomena underlying these di- 
verse data sets requires flexible and accurate models. In this paper, we promote 
Gaussian process d5mamical systems (GPDS) as a rich model class that is appro- 
priate for such analysis. In particular, we present a message passing algorithm 
for approximate inference in GPDSs based on expectation propagation. By posing 
inference as a general message passing problem, we iterate forward-backward 
smoothing. Thus, we obtain more accurate posterior distributions over latent 
structures, resulting in improved predictive performance compared to state-of- 
the-art GPDS smoothers, which are special cases of our general message passing 
algorithm. Hence, we provide a imifying approach within which to contextualize 
message passing in GPDSs. 



1 Introduction 

Analysis and decision making based on data with complex time-varying dynamics 
are of central importance in many areas, such as machine learning, robotics, control 
and economics. The time-series generated in these research areas often have unknown 
dynamical structure, are high-dimensional and noisy, which makes the availability 
of flexible and accurate models highly desirable. We examine Gaussian process dy- 
namical systems (GPDS) as a class of models for exactly such modeling tasks. While 
traditional approaches for time-series modeling and inference, such as the Kalman fil- 
ter [1], often focus on parametric state-space models, the GPDS allows for inference 
over distributions with nonlinear transition and measurement dynamics using Gaus- 
sian processes. 



GPDSs were introduced as Gaussian process dynamical models, a model class that 
is well suited to analyzing high-dimensional human motion in video sequences [19]. 
There are currently two existing approaches for inference in GPDSs: a filtering ap- 
proach based on either a linearization of the GP posterior or particle representa- 
tions of densities [8] and a forward-backward smoothing algorithm using moment- 
matching [6]. These existing approaches are restricted to a single sweep, limiting the 
power of the GPDS model. Iterative algorithms that generalize these filtering and 
smoothing algorithms are desirable as they allow the performance of GPDS inference 
to be improved. 

In this paper, we present algorithms for more accurate inference in GPDSs by devel- 
oping an iterative local message passing framework based on expectation propagation 
(EP). After providing an overview of GPDSs and EP, we make the following contri- 
butions: (1) We develop iterative message passing algorithms for GPDSs based on EP, 
which allows for iterative refinement of the posterior distribution and, hence, better 
inference. (2) We show that the general message passing framework recovers the EP 
updates for a standard dynamical system as a special case and expose the implicit 
modeling assumptions made in standard models. We also show that EP in GPDSs en- 
capsulates all forward-backward smoothers [6] as a special case and transforms them 
into an iterative and more accurate algorithm. (3) We highlight the benefits of our 
EP-based GPDS inference algorithm in the context of a classical tracking problem and 
in the more complex case of inferring low-dimensional posteriors of high-dimensional 
motion capture data. 

2 Gaussian Process Dynamical Systems 

A Gaussian process (GP) is a distribution over functions /, specified by a mean func- 
tion jjLf and a covariance function/ kernel kf [15]. In this paper, we assume that 
the prior mean function equals zero everywhere. Given a set of training inputs 
X = [xi,...,Xn]^ and corresponding training targets y=[yi,... yn]^, yi = f{xi)+ w, 
w ~ M{0,a^), the GP posterior predictive distribution at a test input x^, is Gaus- 
sian distributed M{y* | ;U/(ic*), (7|(cc*)) with mean fJ.f{x^) = kj K^y and variance 
c7y(a;*) = A;** - kj K~^k^ + a^, where fc* = kf{X, a;*), fc** = kf(x^, x^) and K is the 
kernel matrix. In the context of this paper, we assume that the GP models are trained. 

Gaussian process dynamical systems are a general class of discrete-time state space 
models xt = f{xt-i) + wt with wt ~ AA(0, Q) and zt = g{xt) + vt with vt ~ A/'(0, R). 
Here, x G is a latent state which evolves over time, and z e R^, E > D, are 
measurements. We assume i.i.d. additive Gaussian system noise w and measurement 
noise v. The central feature of this model class is that both the measurement function 
g and the transition function / are not explicitly known, but described by posterior 
probability distributions over them. These function distributions are non-parametric 
GPs and we write f ^ QVf and g ~ QVg. Such GPDSs are desirable as they remove 
the restriction of exactly knowing the transition and measurement fimctions by placing 
distributions over them. GPDSs provide a flexible mechanism with which to analyze 
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Figure 1: Factor graph and fully factored graph of a dynamical system. 



the types of rich and complex time-series that are now ubiquitous. 

Due to the nonpar ametric Bayesian properties of the GP, a GPDS is a less restric- 
tive model than dynamical systems based on parametric function approximators (e.g., 
polynomials). Moreover, GPs come with useful analytical properties that can lead 
to good closed-form and efficient approximate inference. Exact Bayesian inference in 
GPDSs is, however, analytically intractable, requiring the use of approximate strategies 
to infer posterior distributions p{X\Z) over latent states X given a set of observations 
Z. Existing approximate inference approaches for filtering and forward-backward 
smoothing and are based on linearization, particle representations and moment match- 
ing [8, 4, 6]. 

A principled incorporation of the posterior GP model uncertainty into a dynam- 
ical system introduces more noise into the inference. In tracking problems where a 
position of an object is not directly observed, this additional source of uncertainty can 
eventually lead to losing track of the latent state. Hence, GPDSs will greatly benefit 
from improved approximate inference methods. In this paper, we develop approx- 
imate message passing using expectation propagation (EP) for GPDSs, i.e., for dy- 
namical systems with posterior GP distributions over the transition and measurement 
functions. We will show that forward-backward smoothing profits from the iterative 
refinement scheme of EP. 



3 Expectation Propagation in GPDS 

EP [10, 11] is an accurate and deterministic approximate inference algorithm, based 
on generalizations of Assumed Density Filtering [3]. EP iteratively refines the pos- 
terior distribution over the latent variables. EP is derived using a factor-graph rep- 
resentation of a graphical model, in which the joint distribution over the latent vari- 
ables p{x) = Hi^il^) is a product of factors /messages ti{x). EP then specifies an 
iterative message passing algorithm in which p{x) is approximated by a distribution 
q{x) = Yl- qi{x), a product of approximate messages qi{x). In EP, q and the messages 
qi are members of the exponential family and q minimizes the KL-divergence KL(p||g). 
EP has been successfully applied to a diverse set of models including sparse regression 
models [17], GP classification [9] and inference in dynamical systems [13, 7, 18]. EP is 
provably robust for log-concave site distributions [17] and invariant under invertible 
variable transformations [16]. In practice, EP has been shown to be more accurate than 
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competing approximate inference methods [9, 17]. 

In the context of this paper, we consider factor graphs that represent dynamical 
systems (see Fig. 1). In these systems, there are three types of messages: forward, 
backward and measurement messages. Furthermore, we assume a fully factored graph 
and are only interested in the marginal posterior distributions p{xi\Z), . . . ,p{xt\Z) 
rather than the full joint distribution p{X\Z) = p{xi, . . . , xt\Z). We consider the case 
that both Xi and Zi are continuous variables and the messages qi are unnormalized 
Gaussians, i.e., qi{x) = Z~^Af{x \ /ij, Sj). Following Fig. 1, we denote messages and 
parameters with respect to forward, backward and measurement messages by symbols 
I>, <1, A, respectively. 

3.1 Expectation Propagation for Approximate Inference 

The EP algorithm is realized practically in two steps: a projection step in which each 
approximate message qi{x) is computed in the context of neighboring factors, and an 
update step in which the posterior distribution q{x) is updated using the new message. 

Factors qi{x) in the graph are updated in the context of all remaining factors. 
First, the cavity distribution q^^^^x) = Ylj^i Iji^) is computed by removing factor qi{x) 
from the joint q{x). Second, an updated factor qf^^{x) is found such that the joint 
qf^^ {x)q\^{x) is a close approximation to ti{x)q\'^{x), where ti is the true factor. The 
close approximation is achieved by finding qf^^ such that the moments of qf^"^ {x)q^'^{x) 
match the ones of ti{x)q\'^{x). In the exponential family, the moments can be com- 
puted by means of the derivatives of the log-partition function (normalizing constant) 
of ti{x)q\''{x) [10, 11, 12]. Alg. 1 summarizes the main steps of Gaussian EP in dynam- 
ical models with a latent-state time series xi,. . . ,xt and corresponding observations 
zi, . . . , Zt- 

3.2 Message Passing in GP Dynamical Systems 

The factor graph and the fully factored factor graph from which we derive our EP- 

based message passing for GPDSs are shown in Fig. 1. For each node/time slice xt we 
compute three messages: a forward, backward and a measurement message denoted 
by qt>{xt), q<tixt) and qA{xt), respectively. To update the posterior q and the message 
qi, i G {A, >, <}, the partition function in Eq. (2) needs to be evaluated. In all cases, 
due to the GP dynamics and measurement models (with nonlinear kernels), the log- 
partition function cannot be evaluated in closed form. We will show, however, that 
for particular approximations, we recover state-of-the-art smoothing algorithms for 
GPDSs. In the GPDS, the partition function we have to evaluate is of the form 

Z-\a) = J p{a\b)q{b)db oc J M{a\ h{b), JlaWib | yu^, J:b)db , (6) 

where g is a cavity distribution and h is a nonlinear function, e.g., the posterior GP 
mean function. Generally, we approximate the partition function deterministically by 
a Gaussian-shaped function in a, i.e., Z~^ M{a \ iJ,^,'Sz). This approximation is 
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Algorithm 1 Gaussian EP for Dynamical Systems 
1: Init: Set all factors qi to AA(0,ooJ); Set q{xi) = p(xi) and marginals q{xt^i) = 

M{0, lO^o/) 
2: repeat 

3: f or t = 1 to T do 

4: for all factors qi{xt), where i = I>, A, <1 do 

5: Compute cavity distribution ^^'(ajt) = q{xt)/qi{xt) = M{xt\ |J^'^,'S\'^) 

with 

sv = i^-^ - Eri)-! , = 5:V(s,- V* - 5:-V.) (i) 

6: Determine moments of ti{xt)q^'^{xt), e.g., via the derivatives of 

\ozZ-\n\\ = log !ti{xt)q\\xt)dxt (2) 

7: Update the posterior q{xt) = M{xt \ Ht, Xlj) and the approximate factor 

qi{xt): 

= /xV + sVyT , ^t = 'E\'-'S\\VlVm-2Vs)'S\' (3) 
V„ := V^\. log Z,ri , V, := V^\i log Zr^ (4) 
gi(a;t) = q{xt)/q\\xt) (5) 

8: end for 

9: end for 

10: until Convergence or maximum number of iterations exceeded 



correct if and only if there is a linear relationship between a and b, i.e., a = Jh, where 
J = ^^^|6=fij is the Jacobian evaluated at 6 = yitf,. We require that J is injective, i.e., 
the dimensionality of h is not smaller than that of a} 

With the implicit assumption of a linear relationship between h and a, we obtain 
the total derivatives d(log Z~^)/d/x^ and d(log Z~^)/dSft of with respect to the 
mean and covariance of the cavity distribution q by applying the chain-rule: 

where = JI4 G R^^-^^^^^ and I4 G E^^-^^-^^^ is an identity tensor. 

When we approximate by a Gaussian-shaped function, we are still free to 
choose a method of computing its mean fj,^ and covariance matrix S2, which also 



^Besides implicit linearization, a number of other approximation strategies exist to compute the mo- 
ments of the product inside the integral in Eq. (2), e.g., sampling-based approximations, using analytically 
tractable lower boimds, or other functional approximations. 
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influences the computation of the Jacobian J. In this paper, we focus on either an 
implicit linearization by moment-matching [14] or an explicit linearization of the pos- 
terior mean function [8], both of which are detailed in Appendix A. Let us now detail 
the derivation of the GPDS messages. 



3.2.1 Measurement Message 

An updated measurement message qi\{xt) is obtained by an application of projection 
to compute the moments of the new posterior marginal q{xt), followed by an update 
using Gaussian division. To apply projection, i.e., computing moments of the product 
of the true factor and the cavity distribution, we compute the partition function 

= j t^ixt)q\''ixt)dxt oc j t^ixt)Mixt\l4\'^}'')dxt, (9) 

t^{Xt) = p{z\Xt) = M{Z I IJig{Xt), ^g{Xt)), (10) 

where tA is the true measurement factor. La Eq. (9), we made it explicit that Z'^^ 
depends on the moments /i}^ and 5]}^ of the cavity distribution. The integral in 
Eq. (9) is intractable due to the nonlinearity of the posterior predictive mean Hyixt) 
and covariance '^gixt) of QVg inn the true measurement factor t/:,{xt) given in Eq. (10): 
Solving Eq. (9) corresponds to a GP prediction with uncertain inputs, which is no 
longer Gaussian [14]. Hence, the integral in Eq. (9) must be approximated. 

We approximate the partition function Z]^^ by a Gaussian Z]^'^ = M{zt \ , Xl^^) 
that depends on the mean fi}^ and covariance of the cavity distribution qt,{xt). 
The parameters of this approximation, a mean fx^^ and a covariance 5]^^ are computed 
analytically: either exactly using moment matching [14], or approximately by lineariz- 
ing the posterior mean function [8]. These parameters are functions of the mean and 
variance of the cavity distribution /i}^ and 5]}^, respectively. The projection step is 
completed by computing the moments via the derivatives of the log-partition func- 
tion: dlogZ^^/d/i)^ and dlogZ^^/dSf , following Eqs. (7)-(8) and (3). The measure- 
ment message is updated by dividing the updated posterior by the cavity distribution 

3.2.2 Backward Message 

To update the backward message q^(xt), for the projection, we require the partition 
function 

Z-\txY,^Y) = j U{xt)q\<{xt)dxt cx j U{xt)M{xt \ fx}^ ,i:)^)dxt , (11) 

t<i{xt) = J p{xt+i\xt)q\:^{xt+i)dxt+i = J M{xt+i\ Hf{xt),'Sfixt)) . (12) 

Here, the true factor i<i(a;t) in Eq. (12) takes into account the coupling between xt and 
xt+i, which was lost in assuming the fully-factored graph in Fig. 1. The computation of 
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the partition function in Eq. (11) remains intractable since the predictive mean iJ.f{xt) 
and covariance 'Ef{xt) of the posterior dynamics GP are nonlinear functions of xt. 
Reordering the integration in Eq. (11) yields 

■^<H/*t'^i^t^) J Q\>ixt+i) j p{xt+i\xt)q\^{xt)dxt+idxt . (13) 

As for the measurement message, we approximate Z^^ by a Gaussian distribution 

using either moment matching [14] or linearization 
of the posterior GP mean fimction [8] for approximating the inner integral in Eq. (13) 
by M{xt+i I /x^^, Xl^^). With this Gaussian approximation of the inner integral, the 
double integral can be solved analytically. The projection step is completed by com- 
puting the derivatives of the log-partition function w.r.t. the moments of the cavity 
distribution q\^{xt) following Eqs. (7)-(8). Using these derivatives, the moments of 
the new posterior q{xt) and the corresponding backward message are updated ac- 
cording to Eq. (3) followed by Gaussian division q<i{xt) = q{xt)/q\^{xt). 

3.2.3 Forward Message 

Similarly for the forward message, the projection step involves computing the partition 
function 

^>HMt+i'^m)=y" t^{xt+i)q^^ {xt+i)dxt+i= j t^{xt+i)N{xt+i\n)X^,^)Xi)^Xt+u 

(14) 

t^{xt+i) = J pixt+i\xt)q\^{xt)dxt = J Af{xt+i \ Hf{xt),'Sf{xt))q\^ixt)dxt , 

(15) 

where again, the true factor t^{xt+i) in Eq. (15) takes into account the coupling be- 
tween xt+i and Xt, see Fig. 1. The moments of the updated (marginal) posterior 
q{xt+i) ~ t|>(£Ct+i)g\|>(a;t+i) are obtained by projection, which is implemented by 

computing the derivatives of the log-partition fimction w.r.t. the moments /x)^^ and 

of the cavity distribution q\^{xt+i). We propose to approximate the true fac- 
tor tt>{xt+i) using a Gaussian approximation q^{xt+i) based on moment matching 
or explicit linearization [14, 8]. We obtain the updated marginal posterior q{xt+i) by 
Gaussian multiplication, i.e., q(xt+i) = qi^(xt+i)q\^{xt+i). 

3.3 Recovering EP Updates for Common Gaussian Smoothers 

The updates of the posterior via the derivatives logZ^^ of our Gaussian approxi- 
mation in the measurement message correspond exactly to the standard EP updates in 
dynamical systems [13], i.e., /x^ = /i}^ + K{zt — fiY) and "Et = — where 
j^xz\A _ (.Qy^xY,z}^] and K = Xl^^^^(I]z^)~^ is the Kalman gain. In a standard 
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single-sweep forward-backward smoother, the cavity distribution corresponds to 
the predictive distribution ]5(a;t|zi;4_i) (time update), and the updated marginal distri- 
bution corresponds to the filter distribution p{xt\zi:t) [1]. In the backward message, the 
posterior updates via the derivatives of the log-partition function log Z^^ correspond 
exactly to common Gaussian EP updates in dynamical systems [13], i.e., fif = fi}^ + 
L{fit+^ - fiY) , and St = + LiSt+i - ^]^)L'^ , where L = cov[x)^ , x^^^ji^E^^^yK 

In a forward-backward smoother (single sweep), the moments /x^^ and of Z^^ 
are the mean and covariance of the time update p{xt+i\zi;t)- The cavity distribution 
q\^{xt) required for the forward message corresponds to the filter distribution p(a;t|zi:t). 
With these direct correspondences, it is straightforward to generalize standard Gaus- 
sian smoothers to iterative EP smoothers. 

Therefore, our iterative message-passing algorithm provides an EP-based general- 
ization and a unifying view on existing approaches for smoothing in (GP)DS, e.g., 
(Extended/Unscented/Cubature) Kalman smoothing and the corresponding GPDS 
smoothers [5, 6]. Computing the messages via the derivatives of approximated log- 
partition functions log as described in Sec. 3.2 recovers not only the standard EP 
updates in dynamical systems [13] and, therefore, also the standard Kalman smooth- 
ing updates. Our message-passing formulation is also general as it comprises all con- 
ceivable (Gaussian) filters /smoothers in (GP)DSs, solely depending on the prediction 
technique used. These predictions define the Gaussian-shaped approximation Z^^ of 
the partition functions and the derivatives of log Z~^ with respect to the moments of 
the cavity distribution, see Eqs. (7)-(8), where the mean and covariance of the predic- 
tion are denoted by /j,^ and Xl^, respectively. 

4 Experimental Results 

We evaluated our proposed EP-based message passing algorithms on three data sets: 
a synthetic data set for illustration purposes, a low-dimensional simulated mechani- 
cal system with control inputs and a high-dimensional motion-capture data set. We 
compared our proposed EP smoothers, EP-GPEKS and EP-GPADS, to state-of-the-art 
forward-backward smoothers in GPDSs, specifically the GPEKS [8] and the GPADS [6]. 
The (EP)-GPEKS uses explicit linearization of the GP posterior mean functions to ap- 
proximate the partition functions, whereas the (EP)-GPADS uses moment matching to 
compute Zf'^. Details of these approximations are provided in the appendix. For the 
synthetic data, we also compared to the Extended Kalman Smoother (EKS) and an EP- 
iterated EKS (EP-EKS). Both GVf and QVg use squared exponential covariance func- 
tions with automatic relevance determination. In all our experiments, we evaluated the 
model using test sequences of measurements Z = [zi, . . . , zt] ■ Whenever available, we 
compared the inferred posterior distribution p{X\Z) of the latent states with the im- 
derlying ground truth using the average Negative Log-Likelihood (NLL^:) and Mean 
Absolute Errors (MAE^). We also report the corresponding average Log-Posterior 
Uncertainty in latent space, which is defined as LPUa; := Z^n=i T~ Z^t=i I'^t I- 
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Table 1: Performance comparison on the synthetic data set. Lower values are better. 





EKS 


EP-EKS 


GPEKS 


EP-GPEKS 


GPADS 


EP-GPADS 


NLL, 


-2.04 ±0.07 


-2.17 ± 0.04 


-1.67 ±0.22 


-1.87 ±0.14 


+ 1.67 ± 0.37 


-1.91 ± 0.10 


MAE^ 


0.03 ± 2.0 X 10"-' 


0.03 ± 2.0 X lO"'' 


0.04 ±4.6 X 10"^ 


0.04 ±4.6 X 10-'^ 


1.79 ±0.21 


0.04 ± 4 X 10-^ 


LPU^ 


-2.66 ±0.03 


-2.62 ±0.03 


-2.51 ±0.04 


-2.50 ±0.04 


1.33 ±0.37 


-2.42 ±0.05 



This error measure evaluates the tightness/ confidence of the posterior. We will show 
that EP inference in GPDS improves typically upon single-sweep forward-backward 
smoothers in all performance measures. 

The computation demand our two presented inference methods is vastly differ- 
ent. High-dimensional approximate inference in the motion capture example using 
moment matching (EP-GPADS) was about two orders of magnitude slower than ap- 
proximate inference based on linearization of the posterior CP mean (EP-GPEKS): For 
updating the posterior and the messages for a single time slice, the EP-GPEKS re- 
quired about 0.5s, the EP-GPADS took about 80s. Hence, numerical stability and 
more coherent posterior inference with the EP-GPADS trade off with computational 
demands. 

4.1 Synthetic Data 

We considered the nonlinear dynamical system xt+i = 4sin(xt) + w , where w ~ 
7^(0,0.12) ^ 4sin(xt) + v, where v ~ M{0,0.l'^). We used p{xi) = Af{0,l) 

as a prior on the initial latent state. We assumed access to the latent state and trained 
the d}mamics and measurement GPs using 30 randomly generated points, resulting 
in a model with a substantial amount of posterior model uncertainty. The length of 
the test trajectory used was T = 20 time steps. Tab. 1 reports the quality of the in- 
ferred posterior distributions of the latent state trajectories using the average NLL^^, 
MAEa; and LPU^ (with standard errors), averaged over 10 independent scenarios. It- 
erated forward-backward smoothing with EP (EP-EKS, EP-GPEKS, EP-GPADS) im- 
proved the smoothing posteriors based on only a single sweep. The standard errors 
in all GP-based methods were relatively high in the NLL^; measure, which is due to 
the remaining model uncertainty. The GPADS had poor performance across all our 
evaluation criteria for two reasons: First, the GPs were trained with few data points, re- 
sulting in posterior distributions with a high degree of uncertainty. Second, predictive 
variances using moment-matching are generally fairly conservative and increased the 
uncertainty even further. This uncertainty caused the GPADS to quickly lose track of 
the period of the state, as shown in Fig. 2(a). By iterating forward-backward smooth- 
ing using EP (EP-GPADS), the posterior was iteratively refined, and the latent state 
could be followed closely as indicated by the small blue error bars in Fig. 2(a) and all 
performance measures in Tab. 1. EP smoothing typically required a small number of 
iterations for the inferred posterior distribution to match the true state. Fig. 2(b). On 
average, EP required less than 10 iterations to converge to a good solution in which 
the average latent-state posterior was essentially identical to the average ground truth, 
see Fig. 2(c). 
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(a) Example trajectory distributions with confi- 
dence bounds. 
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(b) Average NLL^; as a function of the EP itera- 
tion with standard error. 
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(c) Average state trajectories with standard error. 

Figure 2: (a) Posterior state distributions using EP-GPADS (blue) and the GPADS (red). 
Ground truth is shown in black, (b) Average NLL^. per data point in latent space with 
standard errors of the posterior state distributions computed by the GPADS and the 
EP-GPADS as a function of EP iterations, (c) Average trajectories (black) and posterior 
means (blue) with standard errors. 

4.2 Pendulum Tracking 



Figure 3: Performance comparison on the motion capture 
data and pendulum-swing data. Lower values are better. 





Motion Capture 


Pendulum 


Method 


LPU, 


NLL, 


MAE:,. 


LPU, 


GPEKS 


-10.79 ±0.05 


0.29 ±0.30 


0.30 ±0.02 


-2.76 ±0.12 


EP-GPEKS 


-11.01 ±0.05 


-0.24 ±0.33 


0.31 ± 0.02 


-2.77 ± 0.12 


GPADS 


-8.89 ±0.69 


-0.75 ±0.06 


0.29 ± 0.02 


-2.52 ±0.06 


EP-GPADS 


-10.27 ± 0.12 


-0.79 ±0.06 


0.29 ± 0.02 


-2.58 ±0.04 



We considered a pen- 
dulum tracking prob- 
lem in which the state x 
of the system was given 
by the angle (p mea- 
sured from being up- 
right and the angular 
velocity 4>. The pendu- 
lum had a mass of 1 kg and a length of 1 m, and random torques u G [—2, 2] Nm 
were applied for a duration 200 ms (zero-order-hold control). The system noise co- 
variance was set to = diag(0.3^, 0.1^). The state was measured indirectly by two 
bearings sensors with coordinates = (—2,0) and (2:2,2/2) = (—0.5,-0.5), re- 

spectively, according to 2; = [^1,22]''' + , f ~ TV (O, diag(0.1^, 0.05^)) with Zi = 
arctan { ^^s^-x )' ^ ~ trained the GP models using 4 randomly generated 

trajectories of length T = 20 time steps, starting from an initial state distribution 
p{xi) = 7V(0, diag(7r^/16^, 0.5^)) around the upright position. For testing, we gener- 
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ated 12 random trajectories starting from p{xi). 

Fig. 3 summarizes the performance of the various methods. According to the 
LFUx values, iterative methods using EP slightly tightened the posteriors in latent 
space determined by the single-sweep smoothers (GPEKS and GPADS). Generally the 
(EP-)GPADS performed better than the (EP-)GPEKS across all performance measures 
except for the LPUx- This indicates that the (EP-)GPEKS suffered from overconfident 
posteriors than the (EP-)GPADS, which is especially pronounced in the degrading 
NLLa; values with increasing EP iterations and the relatively high standard errors. In 
about 20% the test cases, the inference methods based on explicit linearization of the 
posterior mean function (GPEKS and EP-GPEKS) ran into numerical problems typical 
to linearizations [4], i.e., overconfident posterior distributions resulted in NLL values 
of iboo. We excluded these runs from the results in Fig. 3. The inference algorithms 
based on moment matching (GPADS and EP-GPADS) were numerically stable as their 
predictions are typically more coherent due to conservative approximations of mo- 
ment matching. 

4.3 Motion Capture Data 

We considered motion capture data (from http://mocap.cs.cmu.edu/, subject 64) 
containing 10 trials of golf swings recorded at 120 Hz, which we subsampled to 40 Hz. 
After removing observation dimensions with no variability we were left with obser- 
vations Zi G R^^, which were then whitened as a pre-processing step. For trials 1-7 
(806 data points), we used the GPDM [19] to learn MAP estimates of the latent states 
Xi G M,^. These estimated latent states and their corresponding observations are used 
to train the GP models QVf and QVg. The trials 8-10 were used as test data without 
ground truth labels. The GPDM [19] focuses on learning a GPDS; we are interested in 
good approximate inference in these models. 

Fig. 4(a) shows the latent state posterior distribution of a single test sequence (trial 
10) obtained using the EP-GPEKS posterior distribution. The most significant predic- 
tion errors in observed space occurred in the region corresponding to the yellow /red 
ellipsoids, which is a low-dimensional embedding of the motion when the golf player 
hits the ball, i.e., the periods of high acceleration (poses 3-5). 

Fig. 4(b)-4(c) shows the observed dimensions of the two single time slices with 
the worst and best predictive performance according to predictive performance in 
observed space for the EP-GPEKS method. The horizontal axis shows the observed 
dimensions i = 1 , . . . , 56; the vertical axis shows the values of the corresponding ob- 
servations. The error bars represent the 95% (marginal) confidence intervals of the 
predictive distributions p{z\''\z) = J p{Zf^\xt)p{xt\Z)dxt, i = 1, . . . , 56, and the red 
discs are the actual observations. The worst time slice in Fig. 4(b) is about half-way 
through the motion (slice 77/155), during the high-acceleration phase of the golf swing 
dynamics, see also the yellow/ red ellipsoids between poses 3 and 5 in Fig. 4(a). Except 
for a few dimensions (e.g., 1, 2, 3, 42, 53), the predictions are accurate and confident. 
The overconfidence of some predictions was indicative of a lack of rich training data in 
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(a) Latent space posterior with overlaid golf poses. 




10 20 30 40 50 10 20 30 40 50 

Observation dimension (time slice 77) Observation dimension (time siice 102) 

(b) Observation dimensions of the worst pre- (c) Observation dimensions of the best predicted 

dieted time slices. time slices. 



Figure 4: (a) Latent space posterior distribution (95% confidence ellipsoids) of a test 
trajectory of the golf-swing motion capture data. The further the ellipsoid are sepa- 
rated the faster the movement, (b)-(c): Observation dimensions for two time slices 
(golf data). Red discs are actual observations, the shaded error bars denote the 95% 
confidence intervals of the marginal predictive dimension. 

this area. The best predictive time slice in Fig. 4(c) was obtained when the golf player 
finishes his golf swing just before lowering the golf club (small accelerations). 

Fig. 3 summarizes the results of approximate inference on the golf data set: Iter- 
ating forward-backward smoothing by means of EP improved the inferred posterior 
GPDS distribution. As in the previous experiments, the posterior distribution inferred 
by EP-GPEKS is tighter than the one inferred by EP-GPADS, indicated by the lower 
approximate LPU^: values. NLL^; and MAE^ could not be computed due to the lack of 
low-dimensional ground truth data. 
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5 Conclusion 

We have presented an approximate message passing algorithm based on EP for im- 
proved inference in GP dynamical systems, in order to have iterative and efficient 
message passing algorithms for this important class of models. Our message-passing 
formulation generalizes current inference methods in GPDSs to iterative forward- 
backward smoothing. This generalization allows for improved predictions and has 
been shown to recover existing methods for inference in the wider theory for dynam- 
ical systems as a special case. These special cases are distinguished by the way in 
which the messages are computed. Our new inference approach makes the full power 
of the GPDS model available for the study of complex time-series data. 
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A GP Predictions from Test Input Distributions 

We will now review two approximations to the predictive distribution 



where / ~ QV and Xt-i ~ Af{fi^_i, ^t-i)- 
A.1 Moment Matching 

In the moment-matching approach, we analytically compute the mean fj,^ and the 
covariance of p{xt). Using the law of iterated expectations, we obtain 



where m/ is the posterior mean function of the djmamics GP. For target dimension 
a = I, . . . , D, we obtain 




(16) 



(17) 



V|St-iA-i+/| 



exp ( - ^vJCEt-i + Aa) '^Vi 



) , Vi := {xi - Hi_-i) 



(18) 



for i = 1, . . . , n, where /3„ = K„ 
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Using the law of iterated variances, the entries of for target dimensions a,b = 
1,. . . ,D are 

al = E^,_, [var^[AJ^i_i]] +E;,,,_ JA^] - (19) 
<7^6 = E/,a=,_JA„Aft]-iu>^ a^b, (20) 

respectively, where is known from Eq. (18). The off-diagonal terms u^^ do not 
contain an additional term Ea;j_ Jcovj[Aa, A5|a;t_i]] because of the conditional inde- 
pendence assumption used for GP training: Target dimensions do not covary for a 
given xt-i. 

For the term common to both (7^„ and a^^, we obtain 
E;,,,_, [A„A,] = PlQf3, , Q,j = '^^(-^^'-t-Moj,^.-^) ^1;,Tr^i^^_^^^.) (21) 

with R := (A~^ -|- A^^) -|- / and Zij := A^^i/j -|- A^^uj with Vi taken from Eq. (18). 
Hence, the off-diagonal entries cr^;, of J^t are fully determined by Eqs.(18) and (20). 
From Eq. (19), we see that the diagonal entries (7^„ of "Et contain an additional term 

E,,_, [var/[A„|a;t_i]] = ^ - Tr{K-'Q) + al^ (22) 

with Q given in Eq. (21). This concludes the computation of S^. 

The moment-matching approximation minimizes the KL divergence KL(p||g) be- 
tween the true distribution p and an approximate Gaussian distribution q. This is gen- 
erally a conservative approximation, i.e., q has probability mass where p has mass [2]. 

A,2 Linearizing the GP Mean Function 

An alternative way of approximating the predictive GP distribution for uncertain test 
inputs is to linearize the posterior GP mean function [8]. Given this linearized func- 
tion, we apply standard results for mapping Gaussian distributions through linear 
models. Linearizing the posterior GP mean function yields to a predicted mean that 
corresponds to the posterior GP mean function evaluated at the mean of the input 
distribution, i.e., 

I^t = E/[/a(/*t-i)] = r1/3^ , ra, = al exp ( - ^(x, - A^\xi - (23) 

for i = 1, . . . , n and target dimensions a = 1, . . . , D, where /3(j = K~^y^. The covari- 
ance matrix Ht of the GP prediction is 

= + , V = ^^=Pl^^, (24) 

where Va is given in Eq. (23) and V is the Jacobian evaluated at In Eq. (24), Xl^ 

is a diagonal matrix whose entries are the model imcertainty plus the noise variance 
evaluated at ^^t-l■ This means "model uncertainty" no longer depends on the density 
of the data points. Instead it is assumed constant. 
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Using linearization, the approximation optimality ia the KL sense of the moment 
matching is lost. However, especially in high dimensions, linearization is computa- 
tionally more beneficial. This speedup is largely due to the simplified treatment of 
model imcertainty. 
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